Quantal Andreev Billiards: Density of States Oscillations and 
the Spectrum-Geometry Relationship 
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Andreev billiards are finite, arbitrarily-shaped, normal-state regions, surrounded by superconduc- 
tor. At energies below the superconducting energy gap, single-quasiparticle excitations are confined 
to the normal region and its vicinity, the mechanism for confinement being Andreev reflection. 
Short-wave quantal properties of these excitations, such as the connection between the density of 
states and the geometrical shape of the billiard, are addressed via a multiple scattering approach. 
It is shown that one can, inter alia, "hear" the stationary chords of Andreev billiards. 
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Introduction. The aim of this Letter is to explore cer- 
tain quantal aspects of quasiparticle motion arising in a 
class of mesoscopic structures known as Andreev billiards 
(ABs). By the term Andreev billiard [0 we mean a con- 
nected, normal-state region (N) completely surrounded 
by a conventional superconducting region (S), as sketched 
for the 2D case in Fig. [j]a. The S region is responsible for 
confining quasiparticles that have energies less than the 
superconducting energy gap to the normal region and its 
neighborhood [[y||] . The terminology AB reflects the cen- 
trality of the role played by quasiparticle reflection from 
the surrounding pair-potential ||. Our focus here will 
be on the the density of energy levels of the quasipar- 
ticle states localized near a generically shaped AB, and 
its relationship to the geometrical shape of the AB. The 
main new features that our approach is able to capture 
are the oscillations in the level-density caused by the spa- 
tial confinement of the quasiparticles. This structure is 
inaccessible via conventional quasiclassical methods. 

superconducting 
region (S) 




FIG. 1. (a) Andreev billiard, showing a generic 
retro-reflecting orbit, (b) Orbit corresponding to a station- 
ary chord. The length of this orbit should be contrasted with 
that in (c), which shows a generic "creeping" orbit due to im- 
perfect retro-reflection; electrons (holes) follow full (shaded) 
lines, (d) Generic term in the multiple scattering expansion; 
internal (external) lines represent Green functions G N (G S ). 



Our central results concern the density of states (DOS) 
for billiards of arbitrary shape and dimensionality. They 
include an explicit formula for the coarse DOS, as well as 
a general method for obtaining the (previously inacces- 
sible) oscillations about this coarse DOS, both valid in 
the short-wave limit. This DOS decomposition is feasi- 
ble because, unlike for conventional billiards, the classi- 
cal trajectories of ABs fall into two well-separated classes: 
(i) tracings of stationary chords (see Fig. [lb), and (ii) cer- 
tain extremely long trajectories with many reflections 
(see Fig. |l]c), which contribute to the level density only 
at fine energy resolutions. 

Our strategy for exploring the quantal properties of 
ABs is as follows. First, we express the Green func- 
tion for the appropriate (i.e. Bogoliubov-de Gennes [f|; 
henceforth BdG) single-quasiparticle energy eigenprob- 
lem as an expansion in terms of various scattering pro- 
cesses from the N-S interface. Next, we identify which of 
these scattering processes dominate, by effectively inte- 
grating out processes that involve propagation inside the 
S region, thus arriving at an expansion involving only re- 
flections (i.e. scattering processes that keep the quasipar- 
ticles inside the billiard). The processes associated with 
these reflections can be classified as those that intercon- 
vert electrons and holes (which we refer to as Andreev re- 
flections, and which typically dominate) and those that 
do not (which we refer to as ordinary reflections). We 
then compute the oscillatory part of the DOS via two 
distinct asymptotic schemes. 

The first scheme amounts to an elaboration of that 
adopted by Andreev, and is what is conventionally un- 
derstood when the terms semiclassical or quasiclassical 
are used in the subject of superconductivity. Its physical 
content is perfect retro-reflection (i.e. velocity reversal) of 
quasiparticle excitations from the N-S boundary and per- 
fect e/h (i.e. electron/hole) interconversion (i.e. the ne- 
glect of ordinary reflection processes) . It yields a smooth 
(i.e. low energy-resolution) DOS, as well as singular fea- 
tures that arise from stationary-length chords. However, 
it is incapable of capturing other features in the DOS 
caused by the spatial confinement of quasiparticles. 
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The second scheme incorporates the effect of the im- 
perfectness of retro-reflection which results from differ- 
ences between, say, incident e and reflected h wave vec- 
tors, as well as the effect of ordinary reflection processes. 
It yields the DOS with higher energy-resolution, thus re- 
vealing the oscillations caused by spatial confinement. In 
order to distinguish the effect of imperfect e/h intercon- 
version from higher-order quantum effects, we introduce 
and study a model that features perfect e/h interconver- 
sion but still includes all quantal effects. This model 
is also useful when the pair-potential varies smoothly 
(so that ordinary reflection is even more strongly sup- 
pressed). 

Finally, for the purpose of illustration we examine the 
the case of a two-dimensional circular billiard, and com- 
pare the predictions for the DOS obtained via the various 
asymptotic schemes with those arising from the exact nu- 
merical treatment of the full BdG eigenproblem, as well 
as from the perfect e/h interconverting model. This pro- 
vides a concrete illustration of the implications of wave 
phenomena for the quasiparticle quantum states of ABs. 

Eigenproblem for the Andreev billiard; formulation as a 
boundary integral equation. To address the BdG eigen- 
problem for ABs we focus on the corresponding (2 x 2) 
Green function G, which obeys 
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G(r,r';z) = -I<J(r-r') 



(1) 



where h = — V 2 — k 2 , together with the boundary condi- 
tion that G should vanish in the limit of large |r| . Here, r 
and r' are spatial coordinates, h 2 K 2 /2m is the Fermi en- 
ergy (i.e. k is the Fermi wave vector), h 2 z/2m is the (com- 
plex) energy, and h 2 A(r)/2m is the position-dependent 
superconducting pair potential. The eigcnfunction ex- 
pansion of the Green function leads to the usual repre- 
sentation for the Lorentzian-smoothcd DOS pr(E) of the 
corresponding eigenproblem: 



PT 



* ■* 7T 



r 



(2a) 



ir{E-e n ) 2 + T 2 
- [ lim TrImG(r,r';£ + «r), (2b) 

1" Jr r '^ r 



where Tr denotes a trace over e/h components. 

We assume that the interface between N and S is a ge- 
ometrical surface constituting the boundary of the AB, 
i.e., is perfectly sharp. In other words, A(r) is a constant, 
Ao, outside the billiard and zero inside. Thus, we shall 
not be working self-consistently, but shall benefit from 
being in a position to develop an approach to the quasi- 
particle dynamics that focuses on interface-scattering. 

To construct an expansion for the Green function 
G(r, r'; z) that brings to the fore the geometry of the bil- 
liard (i.e. the spatial shape of the N-S interface) we adopt 
the spirit of the Balian-Bloch approach to the Laplace 
eigenproblem || , and construct a multiple-scattering ex- 



pansion (MSE) in which the Green function is repre- 
sented in terms of the fundamental N or S Green func- 
tions (i.e. those appropriate for homogeneous N or S re- 
gions). Although the physical content of this construc- 
tion is intuitively clear, its development involves lengthy 
technical details which we defer to a forthcoming arti- 
cle The essence of this construction is the deriva- 
tion of a system of integral equations "residing" on the 
N-S interface, the iterative solution of which yields the 
aforementioned MSE for the Green function 0. Within 
this MSE approach, the amplitude for propagating from 
point r in N to r' in N, viz. G(r, r'; z), is expressed as a 
sum of the following processes: (i) the "free" propagation 
amplitude G N (r, r';z); (ii) the amplitude involving a sin- 
gle reflection [i.e. all possible amplitudes for propagating 
from r to a generic interface point a, reflecting at a, and 
then propagating to r': —2 J 9G N (r, a.) cr 3 G N (a, r')]; 
(iii) the amplitude involving two reflections, etc.; (iv) the 
amplitude that traverses the interface twice [i.e. all pos- 
sible amplitudes for propagation from r to the generic 
interface point a, transmission into S, propagation in 
S from a to another generic interface point f3, trans- 
mission into N, and propagation in N from [3 to r': 
-2 2 jf Q 5G N (r, a) <x 3 G s (a, (3) <r 3 5G N (/3, r')]; (v) and so 
on, where a generic term is specified by an ordered se- 
quence of reflections and transmissions (see Fig. [l]d). 
Here, fi,2,3 are the Pauli matrices, and the operators 
d and 5 are defined via 



dG(r,a) 
5G(a,r') 



V r /G(r,r') 
V, G(r,r') 



-ex. J 
-ex J 



(3a) 
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where n a is the normal unit vector pointing into N at ot 
on the N-S interface. 

Semiclassical density of states. So far, our reformulation 
of the BdG eigenproblem has been exact, but many of 
its well-known physical features (such as the dominance 
of charge-interconverting reflection processes) lie hidden 
beneath the formalism. They will, however, emerge when 
we employ either of two distinct semiclassical (i.e. short- 
wave asymptotic) approximation schemes, as we shall 
shortly see. In both schemes, the DOS is calculated via 
Eq. (£§), by using the MSE for G and evaluating the re- 
sulting integrals using the stationary-phase approxima- 
tion, which is appropriate for large kL and small A/k 2 
(where L is the characteristic linear size of the AB). From 
the technical point of view, the difference between these 
schemes lies in the nature of the limits that one assumes 
the parameters to take: 

(A) kL — > oo and A/k 2 — > with LA/ k constant; versus 

(B) kL — > oo with A/k 2 constant. The limit taken de- 
termines which stationary phase points (i.e. classical re- 
flection rules) should be applied. 

In both schemes, however, it is possible to integrate out 
processes involving propagation inside S, to leading order 
in (kL) -1 and A/k 2 . This is done by separating each fac- 
tor of G N and G s in every kernel in the MSE into short- 
ranged pieces and their complements. By doing this we 
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are distinguishing between local processes (i.e. those in 
which all scatterings take place within a boundary re- 
gion of linear size of order so that particles ulti- 
mately leave the boundary region from a point very close 
to where they first reached it), and nonlocal processes 
(i.e. the remaining — or long-range — propagation) . Then, 
we approximate the boundary by the tangent plane at the 
reflection point, and evaluate integrals involving short- 
ranged kernels on this plane. Moreover, contributions 
involving the long-ranged part of G are smaller, by a 
factor of (kL)" 1 , and thus we may neglect them ||. This 
procedure leads to an asymptotic expansion for G, which 
can be used in either of the two semiclassical schemes, 
and which includes only interface reflection (as opposed 
to transmission) and, correspondingly, involves the renor- 
malized Green function G R : 



G~ G f 



dG G 
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G N (r,r')^ fa + (r-r') 
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where tp = cos~ 1 (£'/A), g±(r) = Hq (k±\r\) /4 in two 
dimensions and g±(r) = exp (±ifc±|r|) /47r|r| in three 
dimensions, k± — V k 2 ± E are the e/h wave vectors, 
the integrals are taken over the the interface dV, and, 
e.g., J dv dG N G R = J dv dadG N (r,a)G R (a,r'). Ob- 
serve that the leading term in G R includes only charge- 
interconverting processes; ordinary reflection appears 
only at sub-leading order. In physical terms, the approx- 
imation that we have invoked takes into account the fact 
that an electron wave incident on an N-S interface "leaks" 
into the S side and, consequently, is partially converted 
into a hole and acquires a phase, much as a particle ac- 
quires a phase (i.e. a Maslov index) when reflected by a 
finite single-particle potential. 

We are now in a position to define what we shall call 
the Perfectly Charge-Interconverting Model (PCIM). We 
start with the expansion (||) for G in terms of G R , and 
take the latter to be given by its leading-order form: 
G R w -ie-^erxG^ Then the PCIM is defined via the 
following integral equation for G: 



G = G — 2ie 



'* I dG criG . 

dV 



(5) 



The off-diagonal matrix <j\ ensures that, upon each re- 
flection from the boundary, electrons are fully converted 
into holes (and vice versa). Moreover, this model does 
retain wave propagation effects, as implied by the surface 
integral. 

Let us now focus on semiclassical Scheme A, which 
is, in spirit, the one introduced by Andreev ||. In 
this scheme, excitations undergo perfect retro-reflection 
(i.e. perfect velocity- reversal) , as well as perfect charge- 
interconversion, so that the dynamics is confined to the 



geometrical chords of the AB and, thus, is trivially in- 
tegrable, whatever the shape of the AB Q]. Via this 
scheme, we arrive at the following form for the DOS: 
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Here, the integral is taken over the surface points a 
and /3, and 9 a p denotes the angle between the nor- 
mal at a and the chord leading to /3. This equation 
for pr can be understood as follows: a chord of length 
| a — /3| contributes eigenvalue weight at energies given 
by the well-known semiclassical quantization condition 
En~ 1 \a— f3\ — 2 cos~ 1 (E/A) = 2nn (for n integral). How- 
ever, in order to obtain pr we must sum over all chords 
with the proper weighting, which is accomplished by the 
double integral in over the boundary. The most promi- 
nent features emerging this Scheme A expression for pr 
are singularities, representing the strong bunching of ex- 
act eigenenergies at energies corresponding to stationary- 
length chords. (Such chords have both ends perpendic- 
ular to the billiard boundary.) However, to sum over all 
chords would be superfluous, as the strongest features in 
the DOS can be captured simply from the neighborhoods 
of the stationary chords. Moreover, for finite values of the 
parameters (i.e. kL large but not infinite, and A/ k small 
but non- vanishing) Scheme A produces a locally averaged 
DOS, which becomes numerically accurate only around 
the DOS singularities that it predicts. Thus, it fails to 
capture this DOS oscillations due to the confinement of 
the quasiparticles. The reason for this failure is the fact 
that by summing over all chords one is implicitly assum- 
ing the absence of transverse quantization/confinement. 
To capture such oscillations is the main motivation for 
semiclassical Scheme B. 

In Scheme B we first take into account the imperfect- 
ness in retro-reflection arising from the the previously- 
neglected difference between the wave vectors of inci- 
dent and reflected electrons and holes, whilst neglecting 
all amplitudes involving ordinary reflection. The cor- 
responding classical dynamics is no longer a priori inte- 
grate; on the contrary, it is chaotic for generic shapes [Q. 
In this scheme, the closed periodic orbits fall into two 
classes, quite distinct from one another: one consists of 
multiple tracings of each stationary chord (we refer to 
such chords as As); the other of much longer trajectories 
that "creep" around the billiard boundary (see Fig. |l|c). 
Correspondingly, the DOS is the sum of (i) an average 
term p &v (E), which depends in 3D on the volume (or 
in 2D on the area) of the billiard (i.e. the leading Weyl 
term) ; together with an oscillatory term p osc (E) consist- 
ing of (ii) a finer-resolution term, having a universal line- 
shape that depends solely on the length and endpoint- 
curvatures of the As JoJ , and (iii) very fine resolution 
terms, which depend on the classical dynamics of the bil- 
liard in question: 



Z (E) ~ Re^Z A e 4AA ^ 4 Li d _ 1 _ f (1- 



i(k + -k_)£ A -2itp 
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+ 



jdic orbits 



A po exp iS p 



(6) 



Here, Li n (z) = YlJLi ^ / 3 n is the polylogarithm func- 
tion, <i is the dimensionality of the billiard, w is the 
dimensionality of the degeneracy of the A (e.g. w = 1 
for a circle), Z\ is a slowly- varying real function of en- 
ergy, determining the size of the DOS oscillations, and 
Aa is a measure of the stability of the A, which deter- 
mines whether the "tail" goes towards higher or lower 
energies. For example, an isolated A in 2D would yield 
Aa = sgn(i?i + R2 — £a) — 1 and 
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where R\ and i? 2 are the radii of curvature of the end- 
points of the A [jlO| . The second term in Eq. (0) is the 
contribution from "creeping" orbits (see Fig. [j]c). In it, 
Ap Q is determined by the stability of the orbit, and S, 



po 



is the action corresponding to the orbit. For a typical 
AB, S po > N((k+ - k-)l A - 2<p), where N = (D(n 2 /A) 
and, thus, "creeping" orbits contribute only to the very 
fine details of the DOS. 



bits are the ones that are close to the boundary and, 
for these, consecutive reflections take place very near to 
each other, and thus "see" only the local curvature of the 
boundary. These considerations allow us to perform an 
"adiabatic" approximation to the expansion in Eq. (Q), 
in which we assume that the curvature of the boundary 
varies slowly, relative to the rate at which creeping orbits 
sample the boundary. In Fig. || we compare this adiabatic 
method with the (exact) result obtained by solving the 
full BdG cigenproblem. 



P (E/K) 
30000 



20000 - 



10000 - 



exact 

adiabatic 

w/ normal ref . 




o 

0.066 0.068 0.07 0.072 0.074 0.076 0.078 




FIG. 2. Density of states oscillations for a circular AB: 
kR = 150; A/k 2 = 0.08; smoothing width F/k 2 = 1.1 x 10" 4 . 

For illustration, in Fig. we compare the predictions of 
Schemes A and B with those of the PCIM. The Scheme-A 
result (dashed line) approximates the average behavior of 
the exact DOS for the PCIM (full line). In contrast, the 
Scheme-B result (dotted line) captures the DOS oscilla- 
tions arising from transverse quantization/confinement. 

Thus far in our semiclassical treatment, we have ig- 
nored all amplitudes involving ordinary reflection. For 
non-grazing incidence [i.e. 9 — (tt/2) ~ 1] the amplitude 
for ordinary reflection is very small (~ A/k 2 cos 2 #)). 
However, for orbits that contribute dominantly to the 
oscillatory structure of the DOS, \9 — (tt/2)| <C 1 and, 
therefore, ordinary reflection amplitudes are not negli- 
gible and must be incorporated. This can be done by 
returning to Eq. (^) and re-evaluating the trace formula 
using the full expression for G R (i.e. not just the lead- 
ing, off-diagonal, term). However, these dominating or- 



FIG. 3. Density of states for a circular AB: nR = 150; 
A/k 2 = 0.08; smoothing width F/k 2 = 1.1 x 10" 4 . 

We conclude by emphasizing one particular feature of 
the first term in Eq. (||): this term gives the coarse DOS 
directly, through simple geometrical information in the 
form of the lengths and endpoint-curvatures of the As. 
This feature allows the design of an AB shape that leads 
to a DOS having a predetermined coarse form. More- 
over, as the stationary-chord terms are well separated 
(in time-space) from the creeping orbits, it possible to 
"hear" not only the volume of an Andreev billiard but 
also its stationary chords. 
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of reflections is very large. 

The apparent singularity at E = is an artifact of the 
assumption of imperfectness in retro-reflection; this im- 
perfectness ceases at E — 0. 
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